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ABSTRACT 

A number of previous studies of the fragmentation of self-gravitating protostellar discs 
have involved suites of simulations in which radiative cooling is modeled in terms of a cool- 
ing timescale (fcooi) which is parameterised as a simple multiple (ficooi) of the local dynamical 
timescale. Such studies have delineated the 'fragmentation boundary' in terms of a critical 
value of yScooi (fiait) such that the disc fragments if yScooi < ySa-it- Such an approach however 
begs the question of how in reality a disc could ever be assembled in a state with yScooi < Pmt- 
Here we adopt the more realistic approach of effecting a gradual reduction in yScoob as might 
correspond to changes in thermal regime due to secular changes in the disc density profile. We 
find that the effect of gradually reducing yScooi (on a timescale longer than fcooi) is to stabilise 
the disc against fragmentation, compared with models in which yScooi is reduced rapidly (over 
less than fcooi)- We therefore conclude that the ability of a disc to remain in a self-regulated, 
self-gravitating state (without fragmentation) is partly dependent on the disc's thermal history, 
as well as its current cooling rate. Nevertheless, the eff'ect of a slow reduction in fcooi appears 
only to lower the fragmentation boundary by about a factor two in fcooi and thus only permits 
maximum 'a' values (which parameterise the efficiency of angular momentum transfer in the 
disc) that are about a factor two higher than determined hitherto. Our results therefore do not 
undermine the notion that there is a fundamental upper limit to the heating rate that can be de- 
livered by gravitational instabilities before the disc is subject to fragmentation. An important 
implication of this work, therefore, is that self-gravitating discs can enter into the regime of 
fragmentation via secular evolution and it is not necessary to invoke rapid (impulsive) events 
to trigger fragmentation. 

Key words: accretion, accretion discs - star: formation - gravitation - instabilities - stars: 
formation 



1 INTRODUCTION 



Following the seminal work o fiG anuni there has been con- 

siderable progress in recent years in understanding the behaviour 



siderable progress in recent years in understanding the behavioui 
of self-gravitating accretion dies (see lDurisen et al.l2007 and refer- 



ences therein). A number of simula tions ( Ganimiel2001 



Rice et al] 



I2OO3I : lLodato&Ric3l2004 l2005h have demonstrated that if the 
thermodynamic properties of the disc are evolved according to a 
thermal equation (involving a cooling term parameterised in terms 
of a cooling timescale, fcooi), then the disc may be able to establish 
a self-gravitating, self-regulated state. In this state, the Toomre Q 
parameter: 



CsK 



(1) 



(where c, is the sound speed, k is the epicyclic frequency (equal to 
the angular velocity in a Keplerian disc) and S is the disc surface 
density) hovers at a value somewhat greater than unity over an ex- 
tended region of the disc. Whereas the state Q = 1, corresponds to 
a situation of marginal stability against axisymmetric perturbations, 
in the self-regulated state the disc is instead subject to a variety of 
non-axisymmetric self-gravitating modes whose effect, through the 
action of weak shocks, is to dissipate mechanical energy (i.e. ki- 
netic and potential energy of the accretion flow) as heat. Thermal 
equilibrium is then attained through the balancing of such heating 
by the prescribed radiative cooling: in essence, self-regulation re- 
sults when the amplitude of these modes is able to self-adjust so as 
to maintain thermodynamic equilibrium against the relevant energy 
loss processes. 
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The above studies have all found, however, that such self- 
regulation is only possible in the case that the cooling timescale 
is not too short: stability demands that yScooi = ^cooi/^i ' exceeds a 
critical value wh ich, for discs with adiabatic index of 5/3, is ~ 7 
jRice et alj[20o3) . In the case of more rapid cooling, the disc in- 
stead fragments. 

Such simulations however approach the 'fragmentation 
boundary' in a manner that is unlikely ever to apply to discs in 
reality. In the simulations, the discs are set up without additional 
heating mechanisms and are subject to cooling at some prescribed 
value of P^oo\- Discs with yS^ooi < Pan then fragment on the local 
cooling timescale (i.e. a few times the local dynamical timescale), 
thus begging the question of how such unstable initial conditions 
could ever have been set up in the first place. 

A more likely scenario for disc fragmentation is that the disc 
is instead set up in self-regulated, self-gravitating state and then 
conditions gradually change so that yScooi is lowered. (For example, 
continued infall of material onto a disc or secular re-arrangement of 
material in the disc due to the action of gravitational torques could 
alter the surface density profile of the disc and allow it to enter a 
new cooling regime with lower p^ooi)- It is not however clear that 
the fragmentation boundary would be the same in the case that /Jcool 
is gradually reduced. 

In this paper we conduct a suite of idealised simulations in 
which we explore whether the fragmentation boundary just de- 
pends on the instantaneous value of p^ooi (as has been assumed 
hitherto) of whether the system 'remembers' the history of how 
it evolved to a point of given flcooi- Such a ('toy model') approach , 
is complementary to studies ( Bolev et al fcooel : iMav er et al."2007'; 
Stamatellos et al. 2007, see also the analytical estimates by Rafikov 
20051 : IRafikovl 12007 ) which attempt to achieve ever-increasing 
verisimilitude via the incorporation of more realistic treatments of 
radiative transfer. Here, instead, we make no claims that the sim- 
plified cooling law (for example, the assumption that /3cooi is spa- 
tially uniform) actually corresponds to a situation encountered in 
a real disc, because our aim is to isolate a particular physical ef- 
fect (i.e. the timescale on which the fragmentation boundary is 
approached). The computational expense of 'realistic' simulations 
however prevents their use to study secular effects: even in the case 
of the present 'toy' simulations, it is impracticable to run simula- 
tions over the long timescales on which the E profile changes due to 
gravitational torques or infall. We can nevertheless assess the effect 
of relatively slow changes in jScooi on the fragmentation boundary 
through imposing an ad hoc reduction in the value of /3cooi and can 
apply this insight to the secular evolution of real discs. 

In particular, we want to examine the cause of the fragmenta- 
tion forySeooi < yScrit. that has been found in previous simulations. Is 
this (i) due to the disc's inability to maintain - under any circum- 
stances - a gravitational heating rate that can match the impose d 
high cooling rate? This is the hypothesis o f iLodato & Ric3 j2005l) . 
who identify the minimum value of yScooi with a maximum value of 
the gravitationally induced angular momentum transfer that can be 
delivered by a disc without its fragmenting. They parameterise this 
state of maximal angular momentum transfer in terms of the ratio 
of the r, if) component of the stress tensor to the thermal pressure, 
i.e., by analogy with the equivalent expression for a viscous disc, 
in terms of a m aximum in the well known viscous 'a' parameter 
dShakura & Su nvaev 1973). A critical value of /Jcooi of ~ 7 corre- 
sponds to a maximum or of ~ 0.06. 

Alternatively, (ii) does fragmentation instead reflect the disc's 
inability to set up the required high heating rate on the short 
timescale (fcooi) on which the disc is cooling? If this were the case. 



then with sufficiently gradual approach to the regime of low /3cooi, 
the disc could in principle deliver a value of a that exceeded the 
above limit by a generous margin. 

We can obviously distinguish between these alternatives by 
investigating the case in which pcoo\ is reduced on a timescale t 
that is longer than ?;.ooi, since in this case the disc temperature will 
fall via a sequence of thermal equilibrium states (on timescale r), 
rather than dropping on timescale fcooi- The aim of this investigation 
is thus to see whether the disc is more resistant to fragmentation in 
the regime that t > fcooi • If it is not, then the manner in which the 
disc approaches the fragmentation boundary is unimportant. If, on 
the other hand, it is found that rapid changes in cooling regime are 
required, then it may be necessary to invoke impulsive events (such 
as an external dynamical interaction) to trigger fragmentation. 

In Section 2 we describe the numerical setup, discuss our re- 
sults in Section 3 and in Section 4 we present some conclusions. 

2 NUMERICAL SETUP 
2.1 The SPH code 

Our three-dimensional numerical simulations are carrie d out 
using SP H , a L agrangian hydrodynamic scheme (Ben^ Il990l : 
iMonaghani 1 19921) The g eneral implementation is ve ry similar 
to iL odato & Ricd j2004l) , iLodato & Ric3 j2005l) and iRice et al.l 
(200^. The gas disc is modeled with 250,000 SPH particles 
(500,000 in a run used as a convergence test) and the local fluid 
properties are computed by suitably averaging over the neighbour- 
ing particles. The disc is set in almost Keplerian rotation (allowing 
from slight departures from it to account for the effect of pressure 
forces and of the disc gravitational force) around a central point 
mass onto which gas particles can accrete if they get closer than 
the accretion radius, taken to be equal to 0.5 code units. 

The gas disc can heat up due to pAV work and artificial vis- 
cosity. The ratio of specific heats is y = 5/3. Cooling is here imple- 
mented in a simplified way, i.e. by parameterizing the cooling rate 
in terms of a cooling timescale: 



where u\ is the internal energy of a particle and the cooling 
timescale /cool is assumed to be proportional to the dynamical 
timescale, ?cooi = Pcoo\^ \ where /3cooi is varied according to a time- 
dependent prescription (see Section 2.3 below). 

Artificial viscosity is introduced using the standard SPH for- 
m alism. The actual implementiation is very similar to the one used 
in lRice et alj ( l2005h . that is we set the two relevant numerical pa- 
rameters to ttsPH = 0. 1 and fepn = .2 and we have not included 
here (consisten t with lRice et al.ll200^ the so-called Balsara switch 
l lBalsaralll995h to reduce shear viscosity. 

2.2 Disc setup 

The main physical properties of the d isc at the beginn i ng of the sim - 
ulation are again similar to those o f iLodato & R.ic3 ( |2004 12005[) . 
The disc surface density S is initially proportional to (where R 
is the cylindrical radius), while the temperature is initially propor- 
tional to Rr^^^. Given our simplified form of the cooling function, 
the computations described here are essentially scale free and can 
be rescaled to different disc sizes and masses. For reference, we will 
assume that the unit mass (which is the mass of the central star) is 
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Simulation 


X 


/3hold 


N 


fragmentation 


Fl 


10.5 




250K 


yes 


F2 


10.5 


3 


250K 


yes 


V 


105-10.5 


3 


250K 


yes 


SI 


105 


3 


250K 


no 


Sh 


105 


3 


500K 


yes 


S2 


105 


2.75 


250K 


no 


S3 


105 


2.62 


250K 


yes 


VSl 


314 


3 


250K 


no 


VS2 


314 


2.75 


250K 


yes 



Table 1. Details of the various simulations discussed in this paper. The dif- 
ferent columns indicate: the name of the run, the value of the parameter 
X determining the speed of the reduction of the cooling time, the value of 
yShoM (if any) at which the cooling time was held fixed after reduction, the 
number of particles used in the run N and whether fragmentation did oc- 
cur or not. Simulation V was performed with an initially slow reduction of 
j3 (with X = 105), followed by a fast reduction (with x = 10.5), so that it 
would reach B = 3 with a fast reduction at the same time as simulation S 1 . 



IMq and that the unit radius is lAU. In this units the disc extends 
from /fin = 0.25AU to R„^i = 25AU. The normalization of the sur- 
face density is generally chosen such as to have a total disc mass of 
^disc = 0. IMq, while the temperature normalization is chosen so 
as to have a minimum value of 2 = 2, which is attained at the outer 
edge of the disc. 

Initially, the disc is evolved with constant yScooi = 7.5, this 
value o f jScooi being in th e regime where previous work iGammid 
( l200lh . I Rice et aU jZOOSi) has shown that the disc does not frag- 
ment. T he general features of this initial evolution is described in 
detail in lLodato & Ricd ( l2004h . The disc starts cooling down until 
the vertical scale-length H is reduced such that H/R a Mjisc/M* = 
0.1. At this point the disc becomes Toomre unstable and develops 
a spiral structure that heats up the disc and maintains it close to 
marginal stability. We have evolved the disc with this value of fl^a„\ 
for 7.8 outer disc orbits. At this stage it is close to = 1 over most 
of the disc (i.e. over the radial range R = 3 - 23 A.U. ). 



2.3 Evolution of /3cooi 

After evolution of the disc with fS^ooi = y6cooi(0) = 7.5 for a cooling 
timescale, we effect a linear reduction of /?cooi on a timescale T, i.e. 



j8cooi(0=/Jcooi(0)(l-|:) 



(3) 



where we set T = xfi" '(/?„„,). 

Such a prescription implies that the timescale t(?) on which 
the local instantaneous value of fcooi (i-e. fcooi(^> 0) drops to zero is 



T{t) 



j8cool(0 



R 



(4) 



^cooll /6cool(0) \R, 

We adopt three values of x: x = 10.5 (fast), x = 105 (slow) and 
X = 314 (very slow). In the fast case, t(R, t) ~ tcoo\i.R, t) in the outer 
disc so /(.Qoi is changing faster than the disc can come into thermal 
equilibrium at that value of t^aoi- In the slow case, T(r, t) > tcooiir, t) 
so that the disc is everywhere able to come into thermal equilibrium 
at that This situation is even more amply satisfied in the very 
slow case. 

For each value of x, we run the simulation until a fragment 
forms (at fj^-aoi = jSfias)- I" those cases where the timing of frag- 
mentation suggests that the slow reduction of /J^ooi is acting so as to 
stabilise the disc at lower p^ooi, we test this hypothesis by turning olf 




100 150 

Figure 1. Time evolution of the parameter /3cool in the various models. The 
three solid lines refer to cases where fS^ao[ is first decreased and then held 
at jScool = 3 (that is, simulations F2, SI and VSl). The two dashed lines 
correspond to simulations S3 and VS2. The dotted line is simulation Fl 
and the thick solid line is the higher resolution run (Sh). Finally, the dot- 
dashed line is the simulation with variable rate of change of /Jcool (V). An 
asterisk at the end of the line indicates fragmentation at this time, whereas 
runs without an asterisk are unfragmented at the end of the simulation. 



the reduction in jScooi when it attains a value equal to y6hoid(> Pfmn)- 
We then experiment with values of yShoid in order to find the min- 
imum value of jShoid at which the disc does not fragment over the 
duration of the numerical experiment. Table [T] and Fig. 1 summa- 
rize the main details of the various runs we have performed, where 
the 'F' simulations are the 'fast' ones, the 'S' are the 'slow' ones 
and the 'VS' are the 'very slow' ones (see discussion in Section 3 
below). The simulation named V was performed with an initially 
slow reduction of yS (with x = 105), followed by a fast reduction 
(with X = 10.5), so that it would reach yS = 3 with a fast reduction 
at the same time as simulation Si. This was run as a control run to 
ensure that secular evolution did not alfect our results (see below). 



2.4 Resolution issues 

One important aspect that needs to be taken into account is whether 
the resolution of our simulation is high enough to reproduce frag- 
mentation, when it occurs. Resolution cr iteria for fragm e ntatio n 
with SPH codes have been discussed by i Bate & BurkertI ( Il997l) . 
They obtained that SPH correctly reproduces fragmentation if the 
relevant Jeans mass contains at least 100 SPH particles, that is 
twice the typical number of neig hbours ( A^neigh = 50) within one 
smoothing region. More recently, iNelsonI (^006) has revisited this 
issue focussing on fragmentation in self-gravitating discs and has 
found a slightly more stringent criterion, requiring that the Jeans 
mass is resolved with three times as many particles as required by 
lBate&Burker^ ( ll997l) . In a gravitationally unstable disc, the most 
u nstable wavele ngth is given hy A = 2cJ/GS. The Jeans mass (or, 
as lNelsonll2006l calls it, the "Toomre mass") is then given by: 
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Figure 2. Images of Sh at fragmentation (left) and SI (centre) at the same time. Altliougli Sh has fragmented, only one fragment is seen at large radius and 
aside from the immediate area around the fragment the discs are very similar. This should be contrasted with the profusion of fragments in F2 at the point of 
fragmentation 



M,, = ;rE^2 = ^ = 4;r 



The cumulative disc mass at radius R is given by: 



(5) 



(6) 



since in our setup Z oc ' (see above). We can then rewrite the 
Jeans mass using eq. ([6l(, as: 



Mj 



■■ In- 



(7) 



where we have also used Mdisc = m^Ntat, where A'^t is the to- 
tal number of particles used and m,, is the mass of an individual 
SPH particle. In order to properly resolve fragmentation, we re- 
q uire that Mj > m„N^cso , where N^^^„ = IN^^i^^^ = 100, according 
to iBate & BurkertI jl997l) . or A'reso = SA'ne.gh = 300 according to 
Nelson's more restrictive criterion, and recalling that in our simu- 
lations the mean number of neighbours per particle is 50. We then 
obtain that we have enough resolution at radii R that satisfy: 



R 

^out 



1 



1/3 



(8) 



where q = M^i^^/M-,, and where we have aso used the relationship 
between disc thickness and the parameter Q: 



H 
R 



Q Mi,,,(R) 
2 ' 



(9) 



that can be easilty derived from Eq. l[T]l. Based on equation ([8j we 
can then conclude that for q = 0.1, as used in the present paper, 
fragmentation is well resolved at radii R> 5. Note that, since Meso 
only enters eq. ([8} to th e power of one third, if we had used the more 
restrictive condition of'N elsoril ( 120061 ) ■ we would only increase our 
minimum radius by a factor 1.4. As shown in Fig. 2, whenever we 
observed fragmentation, this occurred outside w 5, so that we can 
be confident that we do resolve the relevant mass and length scales 
for fragmentation. 

A second aspect related to resolution is that we require arti- 
ficial viscosity to play a role only when modeling shocks. In or- 
der to ensure this, we then require that the velocity difference ac- 
cross a smoothing kernel is subsonic, i.e. hQ. < c^, where h is the 
smoothing length. This in turn requires that the smoothing length is 
smaller than the disc thickness H = c^/Sl. We have indeed checked 



that, even at the lower resolution of 250,000 particles, the average 
smoothing length is a fraction » 0.5 of the disc thickness. 



3 RESULTS 

We find that in the fast case {x = 10.5), a fragment forms when 
yScooi = 0.75, i.e. about 8.9 outer disc dynamical times after the 
rapid reduction in y6cooi commenced. Since fragmentation always 
takes about a dynamical timescale to get under way, it follows that, 
as expected, the 'fast' case behaves like the usual case where a fixed 
yScooi is imposed. 

We however see different behaviour in the slow case: here we 
find that when /Jcooi is reduced to, and then held at, jShoid = 2.75 
(the evolution of this simulation is not shown in Fig. 1), the disc 
does not fragment even when the disc is then integrated for a fur- 
ther 44 outer dynamical timescales. Likewise, for the slow case, the 
disc does not fragment when held at yShoU = 3, even after integra- 
tion for 63 outer dynamical timescales at this jScooi value. Q On the 
other hand, when jScooi was instead held at 2.62, it fragmented af- 
ter a further ~ 18 outer dynamical timescales, so it would appear 
that the fragmentation boundary is at around 2.7. This is in strong 
contrast with the value of ~ 7 derived in previous work where a 
fixed y8cooi is imposed. We hesitate to say that we have proved that 
a disc will never fragment when brought to such a low value of 
jScooi value at this slow rate, since our experience shows that where 
one is close to the limit of marginally stable jScooi, fragmentation 
may ensue after long timescales, and that its timing may depend 
on numerical noise that can be affected by resolution. Indeed, we 
found that when we re-ran the x = 105 simulation at higher res- 
olution (A' = 500, 000), and held it at jShoid = 3, it eventually did 
form a fragment at large radius . We however show the disc struc- 
ture in this simulation at the point of fragmentation and contrast it 
with the corresponding situation when the cooling time is rapidly 



' We have tested whether this resistance to fragmentation in the slow case 
is simply because the disc takes longer to reach /3cool = 3 and is therefore of 
lower mass, due to accretion onto the central star However, in the control 
run V (in which the disc attains /3cool = 3 at the same time, but with rapid 
(x = 10.5) reduction in /Jcool between /Jcool = 6 and /3cool = 3), the disc 
fragments promptly. Thus we are satisfied that it is indeed the value of x 
which controls fragmentation. 
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reduced and then held at constant /3cooi = 3 (i.e. model F2). Ev- 
idently, notwithstanding the fact that a fragment does eventually 
form in the former case also, the disc structure is quite different in 
the two cases, with the 'rapid' simulation containing a number of 
regions that are on the point of fragmentation at the moment that 
the first fragment appears. Our interest here is not in defining pre- 
cise boundaries at which fragmentation will or will not occur (since 
the definition of such a boundary is always contingent on the dura- 
tion of the simulation) but in demonstrating that the structure of the 
disc is indeed affected not just by the instantaneous value of /?cooi 
(and hence on the heating rate that has to be delivered through the 
action of the self-gravitating modes) but also on the history of how 
the disc arrived at such a value of yScooi. 

This then raises the possibility (which we discussed in Section 
1) that the lower limit on /Jcooi for self-regulation might represent 
the difficulties that a disc might have in achieving a self-regulated 
state on an appropriately short timescale, rather than a fundamental 
upper limit on the dissipation rate that can be provided by gravita- 
tional modes in the absence of fragmentation. In principle, then, we 
could envisage a situation where the disc might be self-regulated at 
an arbitrarily low value of ^6^001 (i-e- where the non-linear develop- 
ment of the spiral modes delivered an arbitrarily high heating rate 
without the disc fragmenting) provided that the disc approached 
this state sufficiently slowly. Such a conclusion would contradict 
that of Lodato and Rice 2005, who interpreted the fragmentation 
boundary in terms of a (history independent) limit on the maximum 
a value delivered by such instabilities. 

In order to explore this further, we ran the very slow (x = 314) 
simulations so that we could test whether the disc could remain 
self-regulated at a yet smaller value of fj^aoi than for the x = 105 
case. We however found little difference in the results for the x = 
105 and x = 314 case, the lowest values at which yScooi could be 
held being respectively 2.75 and 3.00 for the two cases (for A' = 
250, 000 in both cases). 



4 CONCLUSIONS 

We have found that the rate at which the cooling timescale is 
changed indeed affects the minimum value of jScooi at which the 
disc can exist in a stable, self-regulated state. As expected, this ef- 
fect is only manifest when the the cooling timescale is varied on 
a timescale (r) that is longer than the cooling timescale, since for 
T < ?cooh the temperature always falls on a timescale t^aoi, irre- 
spective of T. We find that when t > t^^oi, the self-regulated state 
is sustainable at cooling times that are about a factor two less than 
those that are possible when a fixed cooling timescale is imposed at 
the outset of the simulation. This implies that (in the slow cooling 
case) the gravitational instabilities are able to deliver about twice 
the heating rate without the disc fragmenting. In terms of the 'vis- 
cous alpha' description of such instabilities (Shakura and Sunyaev 
1973, Gammie 2001, Lodato & Rice 2005), the maximum a de- 
liverable by such a disc is then increased from ~ 0.06 to ~ 0.12. 
It should be noted that such 'local' description of the transport in- 
duced by gravitational instabilities is only possible in the limit in 
which global, wave-like tran sport does not play an important role. 
iLodato & Ricl ( |2004 1200^ , using a cooling prescription similar 
to ours, have shown that this is the case, as long as the total disc 
mass is small (< 0.2M*), which is the case for our simulations. 
iMeiia et al] ( l2005l) , using a constant cooling time, argue that global 
effects might be present, but do not explicitly c alculate such globa l 
torques. On the other hand, recent calculation bv lBolevetai]j2006h 



(see in particular their Fig. 13), which employ more realistic cool- 
ing properties, confirm that in the limit of small disc mass, the 
transport induced by gravitational instabilities is essentially local. 

We have thus found that thermal history can affect the abil- 
ity of the disc to exist in a self-regulated state without fragmen- 
tation but that this affects the location of the stability boundary at 
only the factor two level. The fact that there was negligible change 
in the fragmentation boundary when even slower changes in fcooi 
were employed, demonstrates that thermal history is only part of 
the story. Our results suggest that, however slowly the disc is cooled 
through a sequence of thermal equilibria, there is still a fundamen- 
tal upper limit to the heating that can be provided by gravitational 
instabilities in a non-fragmenting disc. Thus it would appear that 
the initiation of fragmentation in a self-gravitating disc does not 
require that the disc enter the regime of rapid cooling on a short 
timescale. It is thus unnecessary to invoke sudden events (e.g. im- 
pulsive interactions with passing stars, Lodato et al. 2007) to tip 
a previously self-regulated disc into the fragmenting regime. In- 
stead our results suggest that fragmentation can in principle be ap- 
proached via the secular evolution of a self-gravitating disc. 



ACKNOWLEDGEMENTS 

The simulations presented in this work have been performed at 
the UK Astrophysical Fluid Facility (UKAFF). We thank Richard 
Durisen and Ken Rice for interesting discussions 



REFERENCES 

Balsara D. S., 1995, Journal of Computational Physics, 121, 357 

Bate M. R., Burkert A., 1997, MNRAS, 288, 1060 

Benz W., 1990, in Buchler J., ed., The Numerical Modeling of Nonlinear 

Stellai' Pulsations Kluwer, Dordrecht 
Boley A. C, Meji'a A. C, Durisen R. H., Cai K., Pickett M. K., D'Alessio 

R,2006,apj, 651.517 
Durisen R. H., Boss A. P., Mayer L., Nelson A. F., Quinn T., Rice 

W. K. M., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostai's and 

Planets V p. 607 
Gammie C. R, 2001, ApJ, 553, 174 

Lodato G., Meru R, Clarke C. J., Rice W. K. M., 2007, MNRAS. 374. 590 

Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630 

Lodato G., Rice W. K. M., 2005, MNRAS, 358, 1489 

Mayer L., Lufkin G., Quinn T.. Wadsley J., 2007, ApJ. 661, L77 

Mejia A. C, Durisen R. H., Pickett M. K., Cai K., 2005, ApJ, 619, 1098 

Monaghan J. J., 1992, ARA&A, 30, 543 

Nelson A. R, 2006, MNRAS, 373, 1039 

Rafikov R., 2005, ApJ, 621, 69 

Rafikov R. R., 2007, ApJ, 662, 642 

Rice W. K. M., Armitage R J., Bate M. R., Bonnell I. A., 2003, MNRAS, 
338, 227 

Rice W. K. M., Lodato G., Armitage R J., 2005, MNRAS, 364, L56 
Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337 
Stamatellos D., Whitworth A. P, Bisbas T., Goodwin S., 2007, ArXiv 
e-prints, 705 



